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Abstract Correlation functions play an important role 
for the theoretical and experimental characterization 
of many-body systems. In solid-state systems, they 
are usually determined through scattering experiments 
whereas in cold-gases systems, time-of- flight and in- 
situ absorption imaging are the standard observation 
techniques. However, none of these methods allow the 
in-situ detection of spatially resolved correlation func- 
tions at the single-particle level. Here we give a more 
detailed account of recent advances in the detection of 
correlation functions using in-situ fluorescence imaging 
of ultracold bosonic atoms in an optical lattice. This 
method yields single-site and single-atom-resolved im- 
ages of the lattice gas in a single experimental run, thus 
gaining direct access to fluctuations in the many-body 
system. As a consequence, the detection of correlation 
functions between an arbitrary set of lattice sites is 
possible. This enables not only the detection of two-site 
correlation functions but also the evaluation of non-local 
correlations, which originate from an extended region 
of the system and are used for the characterization of 
quantum phases that do not possess (quasi-) long-range 
order in the traditional sense. 



1 Introduction 

The use of ultracold atoms for the study of strongly in- 
teracting many-body systems has undergone remarkable 
development in recent years [1 . Prominent examples in- 
clude the achievement of the strongly interacting regime 
of bosonic and fermionic gases in optical lattices [2- 
|6] and studies of the BEC-BCS crossover by means of 
Feshbach resonances [7 . The success of this approach 
is based on the high degree of control that has been 
achieved over the system parameters. In particular, the 



underlying Hamiltonian is usually known and its para- 
meters, such as the interaction strength or the effective 
mass, can be accurately determined and tuned over a 
large range. Furthermore, ultracold gases experiments 
offer versatile detection techniques, such as in-situ ab- 
sorption, in-situ phase-contrast, and time-of- flight ima- 
ging [4| ISHTO] . The latter has been combined with spec- 
troscopic techniques, such as momentum-resolved radio- 
frequency spectroscopy and momentum-resolved Bragg 
spectroscopy [HI [12], to gain access to the excitation 
spectrum of the many-body system. 
Recent advances in the high-resolution in-situ fluores- 
cence imaging of atoms in optical lattices p!3HT5] have 
pushed detection capabilities to the fundamental level 
of individual atoms. Specifically, this technique allows 
for the single-site-resolved detection of a lattice gas in 
the Bose-Hubbard regime [161 Ull- The Bose-Hubbard 
Hamiltonian is given by 



-J ^ a^jCLi 



u 



^hi{hi - 1) 



'uSi, (1) 



where a] (a^) is the boson creation (annihilation) 
operator on lattice site fii = a\ai is the boson number 
operator, J is the hopping matrix element, U is the 
on-site interaction energy, jii is the local chemical po- 
tential, and the first sum runs over all nearest neighbors. 
The competition of the hopping and interaction process 
gives rise to a quantum phase transition between a 
superfluid and a Mott-insulating phase [TBHTS] . 
The high-resolution fluorescence experiments [131115] 
are able to detect the on-site parity, mod2ni, where 
Ui is the occupation number of site i (see Sec. 2.1). 



The average on-site parity has been used to study 
the superfluid-Mott insulator transition [14 and to 
develop an in-situ temperature measurement in the deep 
Mott-insulating limit [15 . The latter has recently been 
combined with amplitude-modulation spectroscopy. 
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forming a sensitive tool for the study of the many-body 
excitation spectrum, which was used to detect a 'Higgs' 
amplitude mode close to the superfluid-Mott-insulator 
transition |T9l [20]. Furthermore, an antiferromagnetic 
quantum Ising spin chain has been simulated in a tilted 
optical lattice, where nearest-neighbor spin correlations 
map onto the average on-site parity [21 , and an orbital 
excitation blockade could be directly observed [22 . 
The high-resolution technique can also be used to 
control individual atoms and to change Hamiltonian 
parameters on the level of individual lattice sites. In 
particular, it has been experimentally demonstrated 
that the spin of individual atoms in a Mott insulator 
can be addressed [27. This technique has recently been 
employed to study a mobile spin impurity in a strongly 
interacting one-dimensional (Id) system [24] , 
In addition to the average on-site parity, high-resolution 
in-situ fluorescence imaging can be used to evaluate 
correlations between different lattice sites. We used this 
possibility to detect two-site and non-local correlations 
across the superfluid-Mott-insulator transition ^ 
and to image the spreading of correlations after a 
sudden change of the Hamiltonian parameters in the 
Mott-insulating regime [26 . 

Here we give a more detailed account of several 
technical and physical aspects of the detection of 
two-site and non-local correlations with a focus on 
the equilibrium situation in Id [25 . In particular, we 
give a self-contained description of the core experi- 



mental apparatus and sequence (Sec. 2.1). Following 



this, we discuss the observable and the limitation of 



the detection technique in more detail (Sec. 2.2 and 



2.3). Furthermore, we give a short introduction to the 
single-site-resolved detection of Mott insulators in the 
atomic limit [T5] (Sec.[3| as far as it is needed for 
the understanding of the following chapters. After a 



general introduction to correlation functions (Sec. 4.1), 
we demonstrate how two-site correlation functions can 
be used for the detection of correlated particle-hole 



pairs (Sec. 4.2). We give a more detailed introduction to 
two-site parity correlation functions and particle-hole 
pairs than in Ref. [25 and show additional experimental 
data for two-site correlations as a function of the 
in-trap position. We continue with a detailed analysis 
of non-local order in the Bose-Hubbard model (Sec. 
4.3) for which we show a pertubative result that we 
compare to numerical calculations. We finish with 
a short description of the experimental results for 
non-local correlations where our focus is on the role of 
three-site correlations. 



2 Introduction to single-site- and 
single-atom-resolved detection 

2.1 Experimental setup and sequence 

2.1.1 State preparation The starting point of our 
experiments was a two-dimensional (2d) degenerate 
quantum gas consisting of several hundred ^^Rb atoms 
in the hyperfine state 55'i/2 7 F = 1^ mp = —1, which 
was prepared in a single anti-node of an optical lattice in 
the vertical direction (see Fig.jl^). For details concerning 
the preparation procedure, we refer to the supplement- 
ary material of Ref. [25 . 

The vertical lattice was generated by interference 
between an incoming laser beam and its reflection from a 
vacuum window, which has a high-reflectivity coating for 
the laser wavelength Al = 1064 nm, and the lattice spa- 
cing in the vertical direction was aiat = Al/2 = 532 nm. 
The laser frequency was red-detuned to the D2 and Dl 
line, yielding an attractive ac Stark potential [27]. The 
vertical lattice depth Vz was kept constant at approx- 
imately 21 Er, where Er denotes the lattice recoil en- 
ergy Er = / {^maf^^) with m being the atomic mass of 
^^Rb. Due to this tight confinement and an additional 
energy offset from gravity, tunneling of the atoms in the 
vertical direction was negligible for the duration of the 
experiment. Additionally, the vertical confinement was 
much stronger than the harmonic confinement in hori- 
zontal direction that results from the Gaussian beam 
shape. In this pancake-like geometry, the atom cloud 
formed a 2d degenerate Bose gas (see, e.g., Ref. [28] ) . 
Subsequently, the gas was loaded into a 2d optical lat- 
tice formed by two optical lattice axes in the horizontal 
direction, which were created by reflections from mirrors 
outside of the vacuum chamber. The lattice spacing in 
both directions was aiat = Al/2 = 532 nm and the axes 
intersected at 90°, resulting in a simple 2d square lat- 
tice. For the loading, the lattice depths were increased to 
values between bE^ and 23£^r, following s-shaped func- 
tions with durations of 120 ms (Fig.[l]3). These lattice 
ramps were quasi- adiabatic, i.e. slow enough to keep the 
system close to its many-body ground state, as can be 
seen by the low defect density in the system (Sec.|3|. 

2.1.2 Fluorescence imaging We detected the atoms in- 
situ using fluorescence imaging [15, ,29^ . The imaging was 
performed with a high-resolution objective, which was 
placed in front of the vacuum window that is also used as 
a mirror for the vertical lattice axis. The numerical aper- 
ture of the objective was NA = 0.68, yielding a spatial 
resolution of about 700 nm for an imaging wavelength of 
780 nm. 

For the imaging, all lattice depths were increased to typ- 
ically 3000 Er ^ ks ■ 300 fiK on a time scale much faster 
than the many-body dynamics (Fig.[T]3). As a result, the 
original density distribution of the gas was instantan- 
eously frozen. The detailed sequence for the freezing con- 
sisted of a two-step process, where the first step was an 
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Figure 1 Experimental setup and fluorescence imaging, a, Schematic image showing the optical lattice geometry, 
the 2d quantum gas and the high-resolution objective. The optical lattice setup consisted of two horizontal lattice axes (in 
the X and ^-directions) and a vertical lattice axis (in the z-direction). The latter was created by reflection from a vacuum 
window with a coating that reflects the lattice wavelength (1064 nm) and transmits the imaging wavelength (780 nm). The 
2d quantum gas was prepared in a single anti-node of the vertical lattice and is imaged with the high-resolution objective 
using fluorescence detection, b, Horizontal lattice depth during a typical experimental sequence. The loading into the square 
lattice is achieved by quasi-adiabatic s-shaped ramps with a duration of 120 ms. The final lattice depths in x- and ^-directions 
VxjVy were independently controlled. For the detection, the lattice gas is frozen by rapidly increasing the lattice depths to 
approximately 3000 ^r- c. Typical fluorescence image of a dilute thermal cloud as seen on an EMCCD camera. Single atoms 
are visible with a high signal-to- noise ratio. White dots mark the sites of the 2d lattice created by the horizontal lattice axes, 
d. Illustration of the parity-projection mechanism. Figs, a and c are from Ref. [ 15] . 



exponential ramp to ~ SO Er with a duration of 0.2 ms 
that already fully suppressed the dynamics (not shown 
in Fig.[l]3). The second ramp to ^ 3000 Er is included to 
have sufficiently deep lattices for the imaging process. 
To image the frozen gas, we optically pumped the atoms 
from 551/2 7 = 1 to 5Si/2^ F = 2 and shone in ima- 
ging laser light with a red-detuning of ~ 50 MHz to 
the free-space resonance of the 55'i/2, F = 2 to 5P3/2, 
F = 3 transition. The scattered photons were observed 
through the high-resolution objective and an additional 
lens, which focused the light on an EMCCD (electron- 
multiplying charge-coupled device) camera. 
Typical fluorescence images show individual atoms in the 
2d lattice (see Fig.jlJ^). We first discuss a dilute, non- 
degenerate cloud and turn to dense, degenerate clouds 
in Sec.jS] Approximately 5000 photons per atom are de- 
tected during an illumination time of about 1 s, resulting 
in a high signal-to-noise ratio. To achieve this, we had 
to suppress tunneling (or even loss) of the atoms during 
the illumination time. Such tunneling processes can oc- 
cur even in a very deep optical lattice of 3000 E^ if the 
atoms are thermally excited to higher bands as a result of 
the recoil heating by the imaging light. To suppress this 
thermally activated hopping, we laser-cooled the atoms 
by setting up the imaging beams in an optical molasses 
configuration, which resulted in sub-Doppler temperat- 
ures of ~ 30 /iK. 



Note that the imaging process is destructive for the 



many-body state (see Sec. 2.3 for details). To collect stat- 
istics for a given observable, the experiment had to car- 
ried out repeatedly, starting with the preparation of the 
2d degenerate gas. 



2.2 Parity projection and parity operator 

Our images show the parity of the occupation number 
on each lattice site (Fig.[l]i). This is a consequence of 
a pairwise loss process due to light-assisted collisions 
[Ml [151 [30] occurring on a time scale of typically 100 /is, 
which is much faster than the illumination time of about 
1 s. Without additional steps, this loss would occur in the 
beginning of the illumination time, and some of the lost 
atoms would be recaptured due to the molasses cooling, 
leading to an unwanted background signal. To avoid this, 
we trigger light-assisted collisions before repumping and 
before switching on the imaging light with a 50 ms pulse 
of light resonant for the 5Si/2^ F = 2 to 5P3/27 F = 3 
transition. At this stage in the sequence, the atoms are 
in the state 55'i/2, F = 1 and the light is 6.8 GHz red 
detuned for the 55'i/2, F = 1 to 5P3/2, F = 2 transition, 
but efficiently excites atom pairs into molecular states 
leading to a rapid loss. As a result, doubly occupied lat- 
tice sites appear as empty sites in the images. Similarly, 
three atoms on a site are imaged as a single atom as two 
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atoms undergo a pairwise collision. In general, the meas- 
ured occupation number can be written as mod2 in 
terms of the actual occupation number n^. For later use, 
we define an on-site parity operator Si with eigenvalues 
Si that are ±1 for odd (even) parity: 



Si\ni) = Si\ni) 



-\-l\ni) if rii is odd 
— l\ni) if is even. 



(2) 



where \ni) is an on-site Fock state with occupation num- 
ber rij at site i. 



2.3 Description of the observable 

2.3.1 Interpretation as a projective measurement The 
described measurement technique detects more than just 
the average on-site parity but also captures the fiuc- 
tuations and correlations in the system. An arbitrary 
many-body state \^) (not necessarily in the atomic limit) 
at zero temperature can be written as a superposition of 
products of on-site Fock states as 



1^) 



{n^} 



,|ni, ...,nAr), 



(3) 



where |ni,...,nAr) = Hi l^i) o^^^ ^ 

possible configurations of on-site occupation numbers 
{nj = (ni, ...,nAr). 

The freezing of the density distribution and the sub- 
sequent scattering of imaging light can be interpreted as 
a projective measurement. It leads to a projection onto 
a specific state |^)proj = |^i,---,^Ar) with a quantum 
mechanical probability p(ni, tiat) = |Q^ni,...,niv 

2.3.2 Observable without parity projection We first dis- 
cuss the situation as it would be without parity pro- 
jection. The crucial point is that the measurement 
would yield information about all occupation num- 
bers rii of the projected state |l^)proj in a single ex- 
perimental run. In each iteration of the experiment, 
we would observe a new set of occupation numbers 
corresponding to a different \^)proy In this way, we 
could gather an increasing amount of statistics and fi- 
nally reconstruct the full joint probability distribution 
j9(ni, nAr) to observe a specific set of occupation 
numbers (ni,...,nAr) on all lattice sites. This includes 
more information than the average on-site density be- 
cause from the knowledge of p(ni, tiat), one could 
calculate all possible density-density correlation func- 
tions iUieM^i) = E{n,}^K.-.^iv)nieM^i' whcrC 

the product runs over an arbitrarily chosen set of lat- 
tice sites M. 

2.3.3 Observable including parity projection Includ- 
ing the parity-projection mechanism, the most gen- 
eral observable is the joint probability distribution 
Pp(si, Sat) to observe a set of on-site parities 



(5i,...,5Ar) on all lattice sites, where the Si are eigen- 
values of the on-site parity operator as defined in Eq.[2] 
However, the full reconstruction of Pp(5i, sat) is a de- 
manding task since it requires a very large number of ex- 
perimental repetitions. In practice, we evaluated various 
correlation functions using an additional spatial average 
in order to reduce the statistical noise on our data. 

2.3.4 Limitation The detection technique can 
not directly distinguish the pure state \^) = 



E 



|ni, tiat), from a mixed state de- 



scribed by the density operator 



El 

{n^} 



• ,n7v| -P|ni,...,n7v) 5 



(4) 



with 



|ni,...,nAr) 



the projection operator 
|ni, nAr)(ni, nAr|. The reason is that both of 
these states yield the same joint probability distribu- 
tions p(ni, ...,nAr) = |ani,...,n^,P and Pp(si, sat). 
Therefore, the imaging technique, with or without 
parity projection, can not detect off-diagonal elements 
(i.e., coherences) of density operators written in the 
on-site number basis. 



3 Single-site-resolved detection of dense, 
atomic-limit Mott insulators 

3. 1 Theoretical description of the atomic limit 

The atomic limit of the Mott-insulating phase {J /U ~ 0) 
yields a particularly simple description because, neglect- 
ing the tunneling term, the Bose-Hubbard Hamiltonian 
can be written as a sum of on-site terms as 



BH 



u 



hi{hi - 1) - iiihi 



(5) 



The eigenstates of the 
are on-site Fock states 
E{ni) = ^ni{ni - 1) - /HiUi. 
are product states \^)j/u=o 



local Hamiltonian Hi 
\ni) with eigenvalues 
The eigenstates of Hbyl 
— Y\iWi) with eigenener- 



gies E — ^iE{ni). The ground state of the system is 
therefore found by minimizing E(ni) on each lattice 
site. This yields no,i = \iJLi/U^ for the local occupation 
number in the ground state, where \x\ is the ceiling of 
X. The ground state is then 



WW, 



J/U=Q 



(6) 



One of the crucial features of this state is that it shows 
no number fiuctuations, i.e., (n|) — (n^)^ = 0. Number 
fiuctuations only arise through thermal occupation of 
excited states by increasing or decreasing no,*. 
In our system, the Gaussian shape of the lattice beams 
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Figure 2 Imaging Mott insulators in the atomic 

limit, a, Sketch of the phase diagram of the Bose-Hubbard 
model (Eq.[T]) [16 . For small J/U ^ the system tends to be in a 
Mott-insulating phase divided into lobes with constant aver- 
age on-site occupation numbers Ui — {hi) (white areas). For 
J /U larger than a critical value {J/U)c, the system is always 
superfluid (shaded blue area) and only lines with constant r^^ 
exist (dashed lines). Inset: Sketch of the local average density 
of a trapped system in the atomic limit {J/U ^ 0) that can 
be considered as a cut in the Bose-Hubbard phase diagram 
along the /i-direction (grey-shaded arrow), b, Upper panel: 
experimental fluorescence image of a Mott insulator in the 
atomic limit with an occupation number of two in its cen- 
ter. Lower panel: Reconstruction of the atom distribution in 
the lattice. The sites in the central region appear empty due 
to the parity projection mechanism, c, Reconstruction in a 
typical high-density region. Open circles mark the positions 
of single atoms as found by the reconstruction algorithm, 
d, The reconstruction is possible even if the optical resolu- 
tion (about 700 nm) is slightly larger than the lattice spacing 
(ttiat = 532 nm). Five atoms on neighboring lattice sites are 
sketched with their point spread functions (solid line, top 
panel). The sum of the signals (dashed line) shows only a 
broad feature. However, if one of the atoms is missing (lower 
panel) a dip in this feature appears. The reconstruction al- 
gorithm can therefore still distinguish the different occupa- 
tions in the two boxed areas in c. Figs.b,c are from Ref. [15] . 



leads to a harmonic confinement in the center of the 
beams, and the local chemical potential is given by 

1 



(7) 



are 
de- 



where /i is the global chemical potential and uOx^ 
the trapping frequencies in x,?/-directions and x^, y. 
note the x,?/-positions of lattice site i. We do not include 
a confinement term in the z-direction because our sys- 
tem is effectively 2d. The local chemical potential /i^, 
going outwards from the trap center, therefore follows a 
cut ai J/U ~ along the /i axis of the Bose-Hubbard 
phase diagram, starting at the global chemical potential 
/i and going to lower local chemical potentials (Fig.[2|i). 
As a consequence, the Mott lobes appear as concentric 
rings in a single experimental image. 



3.2 Single-site resolved detection in the atomic limit 

Experimentally, we reached the atomic limit by increas- 
ing the lattice depth of both horizontal lattices quasi- 
adiabatically (Sec. |2.1| ) until tunneling between neigh- 
boring sites is completely suppressed. For the data 
presented in Fig.[2j we used a final lattice depth of 
= Vy = 23{l)Er yielding J/U ^ 3 • 10-^ which 
is about a factor of 20 smaller than the critical value 
{J/U)c ^ 0.06 for the superfluid-Mott-insulator trans- 
ition in 2d [3T] . 

After this preparation, we took in-situ fluorescence im- 
ages of the system (Fig.[2]3). The pictures show regions 
with almost constant occupation number and low de- 
fect probability, consistent with the description in terms 
of on-site Fock states. Remaining number fluctuations 
are attributed to thermal excitations of the system (see 
Ref. [15 for a detailed discussion). 

In contrast to the image of a dilute thermal cloud shown 
in Fig.[T]3, single atoms are not easily distinguishable in 
the dense regions of the Mott-insulating shells. However, 
the signal-to-noise ratio and the spatial resolution of our 
imaging procedure are sufficient to reconstruct the atom 
distribution in the lattice using a computer algorithm 
(Figs.[2]3 and d). The output of the algorithm is digital 
information in matrix-form containing the parity of the 
on-site occupation number for each lattice-site J5[ |23] . 



4 Detection of correlation functions 

4-1 Introduction to local and non-local correlation 
functions 

4.1-1 Two-site correlation functions and local order 
parameters The following introduction to local and 
non-local order parameters follows in parts the discus- 
sion in Ref. [32]. Local order parameters are based on 
two-site correlation functions 



(8) 
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where Ak and Bk-\-d are, so far, not specified local oper- 
ators acting at positions k and k-\-d. A system possesses 
long-range order in a two-site correlation function if 



lim C{d) 

d^oo 



|C|Vo. 



(9) 



If this is the case, C is called a local order parameter. 
Typically, one studies the transition from an ordered 
phase, which possesses long-range order, to a disordered 
phase with C = 0. Examples include off-diagonal long- 
range order in a Bose-Einstein condensate defined as 
limd^oo(?^l^fc+d) = I^P 7^ 0, where ^ is the macro- 
scopic wave function and and ^k-\-d are creation and 
annihilation operators. A simple example of magnetic 
order is ferromagnetic long-range order in the Id trans- 
verse field Ising model, defined as 



lim {SlSl^J ^ 0, 

a^oo 



(10) 



where SI. is the z-component of the spin operator. 
In a quasi-ordered phase, e.g., in low-dimensional sys- 
tems, we often find an algebraic decay of C(d) to zero 
and no finite value of C (quasi- long-range order). A dis- 
ordered phase is then characterized by a faster, expo- 
nential decay of C{d). 

4.1.2 Multi-site correlation functions and non-local or- 
der parameters Many phases that do not show 
(quasi-) long-range order in the above sense still show 
non-vanishing long-range correlations in multi-site cor- 
relation functions of the type 



0{l) = (i. 



n 



(11) 



\k<m<k-\-l 



This can be thought of as an extension of the two- 
site correlation function with a string of operators 

n 



k<m<k-\-l 



inserted between the endpoints. We call 
0{l) a non-local correlation function because it gathers 
information about an extended region with length I. A 
system possesses non-local order if 



lim 0(1) = O V 0. 

1^00 

This definition is a generalization of string order 



lim {SI 



n 



^k+h 



^0, 



(12) 



(13) 



\k<ra<k-\-l 



which was introduced in the context of spin-1 chains [33 , 
|34] . This non-local correlation function detects a special 
type of anti-ferromagnetic order, which is hidden for a 
two-site correlation function (Fig.|4|i). Following Refs. 
j32| [35] , we will use the term string order also for the 
more general definition in Eq. ITT] and Eq. 12 



Another example is the disordered, paramagnetic phase 



in the Id transverse field Ising model that shows a non- 
local order 1361 



lim{ n 



(14) 



k<m<k-\-d 



where S^ is the x-component of the spin operator at 
site k. This is in contrast to the ordered, ferromagnetic 
phase of the same model, which shows long-range order 
in a two-site correlation function based on S^ as defined 
in Eq.ITo] 



We will argue in Sec. |4.3| that Id Mott insulators with 
unity filling possess a non-local order given by 



lim ( TT 5^) 7^ 



(15) 



k<m<k-\-l 



with the parity operator Sm defined in Eq.|2] 
This order parameter was introduced to study the trans- 
ition from the Mott insulating phase to a Haldane insu- 
lating phase that can occur if long-range interactions are 
included in the Bose- Hubbard Hamiltonian [33 [38] . The 
latter phase is characterized by an order parameter sim- 



ilar to the one in Eq. 13 In contrast, our focus is on the 



behaviour of the order parameter in Eq. 15 at the trans- 
ition from the Mott insulating to the superfiuid phase 
within the standard Bose-Hubbard model (Eq.[T]). 
Further examples of systems that possess non-local or- 
der are spin- 1/2 ladders [39] and fermionic Mott and 
band insulators [40] . Recent theoretical studies explored 
the connection between non-local order and symmetries 
[35] and the connection to localizable entanglement [411 - 
|44] . A measurement of the string order parameter of 
Eq.[T3] using spin-dependent lattice potentials has been 
proposed in Ref. [45] . 



4.2 Detection of two- site correlation functions and 
particle-hole pairs 

4.2.1 Particle-hole pairs In contrast to the atomic 
limit, Mott insulators at finite J/U > show number 
fiuctuations even at zero temperature and are not de- 
scribed by simple product states as in Eq.[6] For J//7 <C 
(J//7)c, these fiuctuations appear in the form of cor- 
related particle-hole pairs, formed by an extra particle 
and a missing particle on two nearby sites. In the fol- 
lowing, we deal with a situation where the region of in- 
terest lies within the = 1 Mott-insulating lobe and 
the local zero-temperature state, in the atomic limit, is 
described by \^)o,j/u=o — Wi V^i — !)• The emergence 
of quantum-correlated particle-hole pairs can then be 
understood within first-order perturbation theory con- 
sidering the tunneling term as a perturbation. To first 
order, one obtains for the ground state 

|i^)o (X |lZ^)o,j/c/=o + ^ X!«l«il^)o,j/f/=o. (16) 
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The second term yields contributions of the schematic 
form alaj\^)o^j/u=o = V^|l, 1, 0, 2, 1, 1), with a 
sum over all positions and orientations of the particle- 
hole pair formed by the neighboring empty- and doubly- 
occupied sites. The state yields a probability to find 
a particle-hole pair on neighboring sites proportional to 
iJ/Uf. 

Closer to the transition to the superfluid phase, higher- 
order perturbation terms become more important. Intu- 
itively, this leads to a rapid increase of bound particle- 
hole pairs and an extension of their size, eventually res- 
ulting in deconfinement of the pairs at the transition 
point. One might view this proliferation and extension of 
particle-hole pairs as the driving force for the superfluid- 
Mott-insulator transition. However, more complicated 
clusters of particles might also play an important role. 
We would like to stress that \^)o specifies the ground 
state of the system. Particle-hole pairs can therefore be 
regarded as virtual excitations. This should be distin- 
guished from actual excited states of the Mott-insulating 
phase, which are reached by increasing or decreasing the 
atom number by one. 

4-2.2 Two-site parity- correlation functions For the de- 
tection of correlations induced by particle-hole pairs, we 
exploited the fact that they have a finite spatial exten- 
sion. As a consequence, the appearance of a number fluc- 
tuation on a given site leads to an enhanced probability 
to find a fluctuation on a close-by site. This behavior is 
captured in a two-site parity-correlation function [25 , 46^ 



phases [25], which means that {skSk-\-d) ^ {sk){sk-\-d) for 
large d. Therefore, both phases cannot be distinguished 
with Cp{d)^ in contrast to the non-local order parameter 



that we will discuss in Sec. 4.3 We use Cp{d = 1) solely 



(17) 



In contrast to the general definition in Eq.[8j we subtract 
the term {sk){sk-\-d)^ i-e., Cp{d) is a connected correla- 
tion function or a cumulant. In this form, Cp{d) (for 
d ^ 0) vanishes for all many-body states \^) = Yl- {i/ji), 
that are simple products of on-site states because 
such states lead to a factorization {skSk-\-d) = {^k){^k-\-d)' 
Therefore, we expect no signal for the ground state of 
the atomic limit (Eq.|6| and for the ground state of the 
non-interacting limit ([//J ~ 0), which is effectively de- 
scribed by a product of on-site coherent states [4 . Addi- 
tionally, thermally induced number fluctuations do not 
cause a signal in the atomic limit because the density op- 
erator for the finite-temperature ensemble is a product 
of on-site operators [47 . 

In contrast, the perturbation theory result in Eq.[l6 
shows a many-body state that is a superposition of 
product states that have either no particle-hole pair or 
a particle-hole pair at different positions. The nearest- 
neighbor parity-correlation function is then Cp{d = 1) = 
16{J/U)^ + 0[{J/U)^], which is directly proportional to 
the probability to find a particle-hole pair on neighbor- 
ing sites. 

However, Cp{d) cannot be used to formulate an order 
parameter for the superfluid to Mott insulator trans- 
ition because Cp{d) shows a rapid decay with d in both 



to detect quantum fluctuations that are introduced by 
particle-hole pairs in the Mott insulating phase. 

4.2.3 Experimental results We analyzed two-site parity 
correlations in Id systems. For this, we split the 2d 
quantum gas into parallel Id systems by choosing a high 
final lattice depth in ^/-direction of Vy ^ 17 Er effectively 
inhibiting tunneling in this direction (see sequence 
in Fig.[T]3 and illustration in Fig.jS^i). We varied J/U 
within these Id systems by varying the final lattice 
depth in the x-direction from Vx ^ AEr to Vx ~ 17 Er. 
Due to the quasi-adiabatic preparation, we expect the 
system to be close to the equilibrium state for the final 
J/U values. 

Typically, our samples contained 150 — 200 atoms in 
order to avoid Mott insulators of occupation numbers 
rii > 1. For small J/U values, we therefore probed 
the correlations within the = 1 Mott lobe. We first 
evaluated Cp{d = 1) for each nearest-neighbor pair of 
sites in a central region of 9 x 7 lattice sites by an 
ensemble average over 50 — 100 experimental repetitions. 
We then performed a spatial average over this central 
region. The spatially averaged signal can be written 
as - {sk){sk^i))/NR, where Nr is the 

number of nearest-neighbor pairs of sites in the central 
region R. 

First, we recorded the nearest-neighbor correlations 
Cp{d = 1) for different values of J/U along the direction 
of the Id tubes (red circles in Fig.js})). For J/U 0, the 
nearest-neighbor correlations vanish within errorbars, 
compatible with the product-state description of the 
many-body state in this limit. As particle- hole pairs 
emerge with increasing J//7, we observe an increase 
of nearest-neighbor correlations. The appearance of a 
maximum in the signal can be understood from the 
fact that in the non-interacting limit {U/J = 0), Cp{d) 
also has to vanish. The precise location and value 
of this maximum is, however, hard to predict. Note 
that the maximum is reached for a J//7-value that 
is higher than the critical value in mean-field theory 
(J//7)J^p ~ 0.09 but lower than the more precise 
critical value {J/U)\^ ^ 0.3 obtained with numerical 
methods [48", ^49]. For a discussion of the comparison 
with the numerical calculation shown in Fig.jS}), we 
refer to Ref. [25]. Furthermore, we show Cp{d = 1) as a 
function of the distance from the center of the Id tubes 
(see Fig.js]^ and figure caption). 

In addition to the nearest-neighbor parity correlation 
Cp{d = 1), we also evaluated Cp{d) for distances d = 
and d = 2 in our Id systems (see supplementary 
material of Ref. [25]). In contrast to Cp{d = 1), the 
on-site parity fluctuations Cp{d = 0) = {s^) — (sk)'^ 
do not completely vanish in the atomic limit due 
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Figure 3 Nearest-neighbor correlations in Id. a, Top panel: Schematic image of the density distribution for < J/U ^ 
{J/U)c with parallel Id tubes and particle- hole pairs aligned in the direction of the tubes. Bottom panel: Typical experimental 
fluorescence images for J/U = 0.06 (left) and reconstructed on-site parity (right). Particle-hole pairs are emphasized by a yellow 
shading, b, Id nearest-neighbor correlations Cp{d = 1) as a function of J/U averaged over a central region of 9 x 7 lattice 
sites (box in lower right panel of a). The data along the tube direction (x-direction, red circles) shows a positive correlation 
signal, while the signal in orthogonal direction (^/-direction, blue circles) vanishes within error-bars due to the decoupling of 
the Id systems. The curves are first-order perturbation theory (dashed-dotted line), Density-Matrix Renormalization Group 
(DMRG) calculations for a homogeneous system at temperature T = (dashed line) and finite-temperature Matrix Product 
State (MPS) [411 150j calculations including harmonic confinement at T = 0.09 U/ks (solid line). The error bars denote the 
la statistical uncertainty. The light blue shading highlights the superfluid phase, c, Cp{d = 1) as a function of the distance 
6x from the center of the Id tubes (as shown by the arrow in the lower right panel of a). The experimental data (top panel) 
shows a positive correlation signal for intermediate J/U and for < 4 in accordance with the theoretical prediction based on 
MPS simulations (bottom panel). We attribute a negative experimental signal at the edge of the cloud {6x ^ 8) to shot-to-shot 
fluctuations of the cloud size. For the experimental data in c, we averaged the central 7 tubes and additionally averaged over 
data points with the same distance 6x to the left and the right of the center of the Id tubes. The signal in b stems from an 
average over < 4 of the signal in c. Parts of Fig. a, and Fig. b are from Ref. [25] . 



to thermally induced fluctuations, and Cp{d = 0) 
increases monotonically with increasing J/U saturating 
at Cp{d = 0) ^ 1 in the superfluid regime. As a function 
of the distance Cp{d) drops rapidly and the numerical 
calculations predict only a small maximum of 0.01 in 
the next-nearest-neighbor correlation Cp{d = 2) at 
J/U ~ 0.17, which is indiscernible from the statistical 
noise in our measurements f25l. 



4.3 Detection of non-local correlations in the 
Bose-Hubbard model 

4-3.1 Perturbation analysis For the following discus- 
sion of non-local order in the Mott insulating phase, let 
us define a string operator 

dpii)= n (18) 

k<'m<k-\-l 

and the following notation for the string correlator 

olii) = (6p(0)- (19) 

With this, there is string order if 

Ol = Mm Olil) > 0. (20) 

In the atomic limit, at T = 0, the Mott-insulating 
ground state has the form \^)o,j/u=o = Hi 1^* = 



For this state, one finds trivially {Op{l)) = 1 and there- 
fore string order. We will now demonstrate that string 
order is stable with respect to finite tunneling J/U > 0, 
based on the perturbation expansion in Eq.[T6] All states 
that are summed up in the expansion are eigenstates of 
the local parity operator Sm- If no defect is encountered 
at site m, Sm has an eigenvalue 5^ = +1 and, if a defect 
is encountered, we find Sm = —1- Therefore, the states 

are ei- 
= ±1. 



in the expansion on the right hand side of Eq. 16 
genstates of Op{l) with eigenvalues nfe<m<fe+z ' 
The sign of the eigenvalue depends on whether or not 
a particle-hole pair is cut at the ends of the string op- 
erator Op{l). By this, we mean the following: Consider 
a state aja^ |lZ^)o,j/t/=o which has a particle-hole pair at 
positions i and j, where i is the position of the particle 
and j the position of the hole. A particle-hole pair is 
cut if either i or j is within the interval [k^k -\- 1] that is 
sampled by the string operator. If a single pair is cut, the 
eigenvalue of Op{l) is —1, because we find an unpaired 
minus sign. In all other cases, the eigenvalue is +1 be- 
cause particle-hole pairs inside [/c, k -\-l] lead to pairwise 
minus signs that cancel (Fig.|4]3). 

In the expansion of Eq.[T6j we find four different contri- 
butions for which such a cut is possible because there are 
two endpoints of the string operator and two permuta- 
tions of particle and hole. Each of these contributions ap- 
pears with a quantum mechanical amplitude of \f2J /U . 
The total probability to cut a particle-hole pair is there- 
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a Spin 1 chain: hidden anti-ferromagnetic order 

jl^ l //// \ /// 1 1 \ / 
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b Mott insulator: parity order 
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/e<m</e+Z 




J/U 

Figure 4 Illustration of non-local order and numer- 
ical results, a, Hidden anti- ferromagnetic order in a spin-1 
chain. Each spin with ^-component +1 is followed by a spin 
with —1, with an arbitrary number of spins with in between. 
This order is hidden for a two-site correlation function, but 
can be detected using the non-local correlation function in 
Eq.[T3] b, Illustration of parity order in Id Mott insulators. 
For states where all particle-hole pairs lie within the evalu- 
ation region [k, /c + /] , the string operator Y[k<m<k-\-i ^'^n has 
an eigenvalue of +1 because the introduced pairwise minus 
signs cancel (upper panel). Only states which have a particle- 
hole pair extending out of the evaluation region yield — 1 
(lower panel), c, (Pp(/) as a function of J/U calculated with 
DMRG for a homogeneous chain (n = 1, T = 0) of total 
length 216. Lines show Op{l) for selected lengths / (black to 
red colors). The dashed line shows the first-order strong coup- 
ling result Op{l) = 1 — 16(J/?7)^. Inset: Extrapolated value 
Op = lim^^oo Op{l) together with a fit (black line) of the 
form (9^ oc exp ( - A [{J/U)c - {J/U)]~^^^), which is char- 
acteristic for a transition of Berezinskii-Kosterlitz-Thouless 
type. Fig. c without the strong coupling result was first pub- 
lished in Ref. [25,. 



fore Pc = A{^/2J/Uf = ^{J /Uf and the total probabil- 
ity not to cut a particle-hole pair is Pnc — 1 — S{J /UY . 
Based on the argument of the previous paragraph, we 
can write 

Ol(l)=Pnc-Pc = l-l^{J/U)\ (21) 

This result is valid if / > 2 and for J/U <^ {J/U)c. Con- 
sequently, we find \imi^oo{Ol{l)) = 1 - 16{J/U)'^ > 0. 
Therefore, Mott insulators in Id at small J/U possess 
string order. 

The previous analysis is expected to hold in a similar 
way in the whole Mott phase if higher order terms in 
the perturbation series are included. The crucial point 
is that for a given order (J/U)^ in the Mott-insulating 
phase, defects appear in finite-sized clusters with a typ- 
ical extension /c- For low orders of n, we expect Ic ^ n. 
Within such clusters, the number of lattice sites with 
eigenvalue = — 1 is always even because the applic- 
ation of the perturbation operator aja^, with z,j be- 
ing neighboring sites, can only change this number by 
±2. For example, the application of aja^ on a state 
|1, 1, 1, 0, 2, 1.., 1, 1) can create another particle- hole 
pair, destroy the particle-hole pair, or lead to states 
|1,1,...,0,1,2,...,1,1) or |1,1,..., 0,3,0,..., 1,1). In all 
cases, we find an even number of sites with s^n = — 1- If 
such a finite-sized cluster lies completely within the eval- 
uation region [/c, + it contributes only with an overall 
plus sign since the contribution of all sites with Sm = 
cancels pairwise. The string-correlation function 0^(1) 
then takes on a constant value if I Ic because clusters 
can only be cut at the ends of the evaluation region and 
the probability for such a cut stays constant as a func- 
tion of the length 

4-3.2 Numerical analysis This intuitive view is suppor- 
ted by a numerical analysis using DMRG (Fig.[4]3, solid 
lines). The result from the strong coupling expansion 
Ol ^ l-16(J/[/)2 (dashed line) agrees with the DMRG 
data for small J/U values. The inset of Fig.[4]3 shows the 
extrapolated values Op — lim^^ooOplO? which we cal- 
culated using finite size scaling (see supplementary ma- 
terial of Ref. [25] for details). We observe that (9^ > in 
the Mott phase and that O^ vanishes at J/U ~ 0.3, 
which is compatible with the known numerical value 
of {J/U)c ~ 0.3 for the critical coupling strength for 
the Id superfluid-Mott-insulator transition [48, 49 . This 
analysis shows that Id bosonic Mott insulators possess 
string order for all coupling strengths J/U < {J/U)c' 
Without going into details, we note that an analysis us- 
ing a Bosonization treatment yields an algebraic decay 
to zero of Op{l) with in the superfluid phase, and 
an approximately exponential decay to a constant non- 
zero value, in the Mott-insulating phase ^J. Further- 
more, a mapping to a 2d interface model shows that 
Ol grows as O^ (x exp { - A [{J/Uy/ - {J/U)]~^^^ ) 
at the transition point and our numerical data is 
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J/U 

Figure 5 Experimental results for non-local correla- 
tions, a, Experimental values of Op{l) for lengths < / < 8 
b, Our experimental values for three site cumulant (siS2S3)c 
(green circles) show the existence of three-site correlations 
in our system. The curves are DMRG calculations for a 
homogeneous system at T = (dashed line) and finite- 
temperature MPS calculations including harmonic confine- 
ment at T = 0.09 [///cs (solid line). Figs. a,b are from 
Ref. EH. 



site cumulant {siS2S^)c^ defined as [54] 

{siS2h)c ={siS2h) - {SI){S2){S2,) 

- C^(l, 2){h) - Cp(2, 3)(^i) - C^(l, 3)(^2), 

(22) 

with two-site correlation functions Cp{i^j) = (siSj) — 
{si){sj). The three-site cumulant {siS2Ss)c is a measure 
of the correlations between all three sites as it vanishes 
if one of the sites is not correlated with the others [ 54] . 
If {siS2Ss)c vanishes, Eq. 22 can be written as {siS2Ss) = 
(5i)(s2)(53)+Ci,2(s3) + C2,3(5i) + C'i,3(52). In this casc, 
{siS2Ss) does not contain more information than two-site 
and on-site terms. In contrast, a non-zero {siS2Ss)c sig- 
nals that (S1S2S3) contains additional information. Our 
experimental values for the three-site cumulant {siS2Ss)c 
(FigjsjD) show a significant signal for J/U < 0.20, in 
quantitative agreement with the MPS calculation. 
In summary, for J/U ~ 0, the string-correlation sig- 
nal is purely dominated by on-site terms while for 
^ J/U ^0.1, also two-site correlations (quantified by 
the cumulant in Eq. 17) start to play a role. In the range 
0.1 < J/U < 0.2, three-site correlations (quantified by 



the cumulant in Eq. 22 ) build up and also contribute 



to the string-correlation signal. In general, one expects 
higher-order correlations to play a more significant role 
in the vicinity of the phase transition. A study of cu- 
mulants including more than three sites is a matter of 
future investigations. 



5 Discussion and outlook 



compatible with this behavior (inset of Fig.jlj^). The 
Id superfluid-Mott-insulator transition is of Berezinskii- 
Kosterlitz-Thouless type, and the latter scaling is the 
same scaling as for the opening of the Mott-insulating 
gap at the transition point [53] . 

4-3.3 Experimental results Our experimentally ob- 
tained values of Op{l) for string lengths / < 8 are shown 
in Fig. [5^. We observe a stronger decay of 0^(1) with / 
compared to the T = case because at finite temperat- 
ure, thermal fluctuations lead to minus signs at random 
positions and reduce the average value of 0^(1). Des- 
pite this, we still see a strong growth of Op{l) once the 
transition from the superfluid to the Mott insulator is 
crossed, similar to the behavior in Fig.[4|3. For a more 
detailed comparison with numerical results, we refer to 
Ref. [25 . 

Here we would like to stress that this signal cannot be 
explained with pure two-site correlations. To illustrate 
this, let us consider the simplest case of a string-type 
correlator including three sites (5152S3), where si,52 and 
ss refer to the parity of three neighboring sites on a Id 
chain. Three-site correlations are captured by a three- 



We presented experimental and theoretical details for 
the single-atom and single-site-resolved detection of cor- 
relation functions in a strongly interacting Id lattice sys- 
tem [25] . More specifically, we studied two-site and non- 
local parity correlations across the Id superfluid-Mott- 
insulator transition. 

In Ref. [25] , we also showed the detection of two-site cor- 
relation functions across the 2d transition. It is an inter- 
esting question whether non-local correlation functions 
can also be used to detect the phase transition in 2d 
using a function Op{A) = (YiieA^'i) ^^^^ evaluates the 
product of on-site parities on an area A. In contrast to 
the Id situation, the boundary dA of the evaluation area 
grows when increasing the area. One can show that even 
for first order in J//7, this leads to an exponential decay 
log [(9^ (A)] (X —dA in the Mott-insulating phase. How- 
ever, the superfluid phase is expected to show a faster 
decay with log[0^(A)] oc —dAlog{dA), and the phase 
transition should be detectable by these different scal- 
ings [52 . 

Furthermore, the detection of single-atom and single- 
site-resolved correlation functions is a valuable tool for 
the study of out-of-equilibrium behavior. In Ref. [26 , we 
could detect the spreading of correlations after a quench 
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in the Id Mott insulating regime. It would be interest- 
ing to extend these studies to 2d systems, where the 
numerical calculation of out-of-equilibrium dynamics is 
demanding. 

Another future direction is the single-site and single- 
atom-resolved detection of ultracold fermions in the 
Fermi-Hubbard regime. The measurement of correla- 
tions could then be used for a direct detection of anti- 
ferromagnetic order at low temperatures. 
Additionally, in a recent experiment, we used the same 
technique to detect correlations between Rydberg atoms 
forming spatially ordered structures [55 , showing that 
the general approach is not restricted to Hubbard- 
type physics. Furthermore, using off-resonant Rydberg- 
dressing [5S^ , an extended Bose-Hubbard model with 
longer-range interactions that features a Haldane insu- 
lating phase [ST, 'SS could be realized. This phase is ex- 
pected to show a non-local order similar to the hidden 



anti- ferromagnetic order in Eq. 13 



The possibility to measure higher-order correlation func- 
tions offers a precise way of characterizing different 
quantum phases. In the long term, however, one wishes 
to perform measurements that are not restricted to 
density-type correlation functions in order to detect 
coherences in the many-body state (see discussion in 



Sec. 2.3). Therefore, a major experimental step would 



be the single-site-resolved detection of the single-particle 
density matrix {a\aj). Based on this, a long-term goal 
would be the direct measurement of correlators involving 
higher orders of creation and annihilation operators that 
do not resemble density correlations. One particularly in- 
teresting application is the investigation of entanglement 
in many-body systems [SH]. Indeed, recent proposals sug- 
gested that single-site- and single-atom-resolved imaging 
could be used for the detection of entanglement in bo- 
sonic [53 and fermionic [60 systems. Possible applica- 
tions range from the detection of entanglement growth 
after quantum quenches to experimental investigations 
of area laws. 
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